function [acc, m, n] = heateqfinitediff(x0, xn, dx, t0, tn, dt)

% [acc] = heateqfinitediff(x0, xn, dx, t0, tn, dt)
%
%   Solves the heat equation diffusion problem in 1d using finite
%   differences on the domain x = [x0, xn] from t = [t0, tn].
%
%               (-4*(x-5)^2)
%   F(x, 0) = e^
%
%   fx0 = initial temperature function
%   dx = spatial step
%   dt = time step

syms fx0x;
fx0 = inline('exp(-4*(x-5)^2)');

m = 1 + (tn - t0) / dt; % (j) y-axis
n = 1 + (xn - x0) / dx; % (i) x-axis

acc = zeros(m+1, n+1);

for i = 1:n
    acc(1, i) = feval(fx0, dx*(i-1));
    acc(m+1, i) = dx*(i-1);
end

for j = 1:m
    acc(j, 1) = 0;
    acc(j, n) = 0;
    acc(j, n+1) = dt*(j-1);
end

for j = 2:m
    for i = 2:(n-1)
        acc(j, i) = acc(j-1, i) + (dt/dx^2)*( acc(j-1, i+1) - 2*acc(j-1, i) + acc(j-1, i-1));
    end
end